Laminar flow of two miscible fluids in a simple network 
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When a fluid comprised of multiple phases or constituents flows through a network, non-linear 
phenomena such as multiple stable equilibrium states and spontaneous oscillations can occur. Such 
behavior has been observed or predicted in a number of networks including the flow of blood through 
the microcirculation, the flow of picoliter droplets through microfluidic devices, the flow of magma 
through lava tubes, and two-phase flow in refrigeration systems. While the existence of non-linear 
phenomena in a network with many inter-connections containing fluids with complex rheology may 
seem unsurprising, this paper demonstrates that even simple networks containing Newtonian fluids 
in laminar flow can demonstrate multiple equilibria. 

The paper describes a theoretical and experimental investigation of the laminar flow of two misci- 
ble Newtonian fluids of different density and viscosity through a simple network. The fluids stratify 
due to gravity and remain as nearly distinct phases with some mixing occurring only by diffusion. 
This fluid system has the advantage that it is easily controlled and modeled, yet contains the key in- 
gredients for network non-linearities. Experiments and 3D simulations are first used to explore how 
phases distribute at a single T-junction. Once the phase separation at a single junction is known, 
a network model is developed which predicts multiple equilibria in the simplest of networks. The 
existence of multiple stable equilibria is confirmed experimentally and a criteria for their existence 
is developed. The network results are generic and could be applied to or found in different physical 
systems. 
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I. INTRODUCTION 



When a fluid comprised of multiple phases or constituents flows through a connected fluidic network, it has been 
observed in different applications that the phase distribution may be exhibit unsteady or non-unique flow for fixed 
inlet conditions. Such heterogeneous distribution of phase within network flows has been studied at a variety of 
scales. At the micro-scale, the flow of droplets or bubbles through microfluidic networks can demonstrate bista- 
bilty and spontaneous oscillations [H-IH- These network non-linearities have been exploited by researchers who have 
demonstrated microfluidic memory, logic,_and control devices On the macro-scale, models of magma flow with 

either temperature-dependent viscosity @ or volatile-dependent viscosity 0] have shown the existence of multiple 
equilibrium solutions on the pressure-flow curve which can lead to spontaneous oscillations. 

Another network that can exhibit complex behavior is micro- vascular blood flow. Nobel prize winner August Krogh 
noted the heterogeneity of blood flow in the webbed feet of frogs in the early 1920's Q. In the Anatomy and Physiology 
of Capillaries he wrote Q 

In single capillaries the flow may become retarded or accelerated from no visible cause; in capillary anas- 
tomoses the direction of flow may change from time to time. 

Numerous researchers have confirmed these observations over the years. De Backer et al. reported the first visualiza- 
tion of micro-vascular flow in critically ill patients and observed alterations in red blood cell distribution in patients 
with sepsis [l(| • Their videos (supplementary materials) capture the heterogeneity of blood flow that Krogh referred 
to; spontaneous changes in flow speed and spontaneous change in the direction of flow in loops. 

The distribution of red blood cells in micro- vascular blood flow arc often interpreted as evidence of biological control. 
If the flow in a vessel increases, it is assumed that the diameter of the vessel responds in order to auto-regulate the 
flow. This vasomotion is assumed to be the cause for oscillations in the micro-circulation [TTj] . While the importance 
of vasomotion cannot be denied, there is significant evidence that fluctuations in cell distributions in micro- vascular 
networks can be due to inherent instabilities [HI, [l3j]. In 2007 Geddes et al. demonstrated that two important 
rheological effects are necessary for the existence of multiple equilibria and spontaneous oscillations in microvascular 
networks — the Fahraeus-Lindqvist effect, which governs viscosity of blood flow in a single vessel, and the plasma 
skimming effect, which describes the separation of red blood cells at diverging nodes [T3 |. 




The rheology of blood has been well studied, with the first comprehensive measurements conducted by Fahraeus 
and Lindqvist in 1931 [l5|. While blood rheology has many complicated details, for the existence of multiple equilibria 
the only required ingredient is that viscosity is a non- linear function of red blood cell concentration (hematocrit) flij . 
While the role of non-linear viscosity on network flows was perhaps first discussed in the context of blood flow, non- 
linear viscosity is a common feature of essentially all fluid systems involving multiple phases. In a recent experiment, 
we (JBG and BDS) demonstrated that multiple equilibria exist in a simple network involving two miscible and well 
mixed Newtonian fluids of different properties due to non-linear viscosity [lfj . These experiments used sucrose solution 
and water, demonstrating that for existence of multiple equilibrium in the pressure-flow curves, there is no need for 
complex rheology. 

Plasma skimming is another source of heterogeneity in microvascular networks [bfj . Krogh introduced the term 
plasma skimming in 1921 in order to explain the disproportionate distribution of red blood cells observed at single 
vessel bifurcations in vivo @. In the absence of plasma skimming, the hematocrit in two downstream vessels of a 
bifurcation would equal that of the feed vessel. Numerous authors have demonstrated plasma skimming in vitro and 
in vivo and many attempts to measure the plasma skimming function have been made p7M22l | . These measurements 
have led to empirical relationships which are a critical component of the theory and simulation of microvascular flow. 

Plasma skimming (or more generally phase separation) exists in numerous other fluid systems. In two fluid systems, 
it is commonly observed that the phase fraction after a diverging bifurcation is different than that in the feed. Probably 
the most widely studied example of such phase maldistribution is gas-liquid two phase flow. This system has important 
technological applications in power and process industries and oil production. In many process applications phase 
maldistribution can have detrimental consequences for downstream equipment [23], while in some cases the phenomena 
is exploited to build simple phase separators [24j . Extensive experimental work on gas liquid flow has been conducted 
over the past 50 years and these studies have been extensively reviewed [23], [2^, [2(| . All this work has shown that 
significant separation can occur which is a function of the inlet volume fractions, the geometry of the bifurcation, the 
fluid properties, and the two-phase flow regime. Simple models can capture some of the experimental features of this 
system |27| . 

In applications for the process and petroleum industry, phase separation in liquid-liquid flows are less well-studied 
though several recent papers have emerged (28|-|30|. These studies have typically been conducted with immiscible 
fluids and at high Reynolds number where the flows are turbulent and often the state of the incoming two-phase 
flow to the bifurcation is critical to the phase separation. Similar phase separation effects occur in systems with 
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liquid- vapor flows with evaporation or condensation. The impact of phase maldistribution in two-phase flow has been 
shown to impact network flows in refrigeration systems [3l| and solar power systems [32| . 

In this work we continue a systematic theoretical and experimental investigation of heterogeneity in simple network 
flows involving ordinary fluids [HI ]. Our aim is to begin to understand the underlying mechanisms that govern the 
phase distribution within networks involving fluids with more than one constituent. We focus on laminar viscous flow 
of two miscible fluids with different viscosity and density. The fluids stratify in the system due to gravity and remain 
as nearly distinct phases with some mixing occurring only by molecular diffusion. This fluid system has the feature 
that it is easily controlled and modeled, yet has the key ingredients for complex network flows. 

In this paper we first use experiments and 3D Navier-Stokes simulations to explore how the two fluids distribute 
at a single T-junction. While such phase separation functions have been extensively measured for gas-liquid flows, 
to the best of our knowledge they have not been measured in this miscible laminar two-fluid flow. We measure the 
phase separation at a single junction, find excellent agreement with 3D Navier-Stokes simulations, and propose a 
simple parametric phase separation function. Once the behavior of a single junction is characterized, we construct a 
simple network model. We find that the phase separation which occurs at a T-junction can lead to multiple stable 
equilibrium in even the simplest of networks. Our experiments confirm these predictions. 



II. PHASE SEPARATION FUNCTIONS 

The first step toward predicting the distribution of phases in a network is to understand the phase separation at a 
single bifurcation. In microvascular flow, these "plasma skimming" functions have received much attention and are a 
crucial component of network models. In gas-liquid flow, these phase separation functions have also been extensively 
measured. In order to make predictions on networks with our fluid system, these phase separation functions must 
first be measured. Our focus is on laminar stratified flow of Newtonian fluids with different viscosity and density. To 
the best of our knowledge these phase separation functions have not been measured for this system. The functions 
are complicated because the depend upon many parameters in the system. In the first part of this paper, we seek to 
explore the phase separation at a single T-junction both experimentally and computationally. 

A. Experimental system 

A top view schematic of our experimental setup is shown in Figure [TJ in this figure gravity points into the page. 
Two syringe pumps (New Era Pump Systems NE-300) supply our source fluids at a controlled and steady flow rate. 
Inlet pump I contains water (which will be denoted as fluid I) and inlet pump 2 contains a controlled aqueous sucrose 
solution (fluid 2). The mass fraction of sucrose (Sigma Aldrich) in fluid 2 is precisely measured with an analytical 
balance when the solution is prepared. The viscosity and density of the inlet sucrose solution are taken from the CRC 
Handbook from the known mass fraction (33|. Circular tubing (1.6 mm inner diameter) from the two inlet pumps 
meet at the inlet junction where the density difference of the two fluids is sufficient to create a strongly stratified 
flow. The dense sucrose solution is observed to sit on the lower half of the tube and the water sits on top. Food 
coloring is added to both fluids for basic visualization. The inlet tube then approaches the test T-junction which is 
held level with respect to gravity (see Figure [T]) . In the schematic we adopt the nomenclature for the two outlets as 
"branch" and "run" ; terms we will refer to throughout the paper. The flow rate of one of the outlets is controlled 
with the outlet syringe pump set in withdraw mode (Harvard Apparatus, ml pump module). The other outlet is held 
at atmospheric pressure and its flow rate is known from the difference of the pump flow rates. The fluid leaving the 
system from the open outlet is collected. The collected fluid can then be measured with a Brixometer (Ataga, PAL-1) 
to determine the total amount of sucrose in the outlet fluid. The volume fraction of fluid 2 in the outlet, can then 
be determined from this measurement and knowledge of fluid contents of the inlet pumps. The volume fraction is 
defined as $ = V2HV1 + V2), where V is the volume and the subscript refers to which fluid. 

The test T-junction has a circular cross section of 1.1 mm inner diameter. The entrance and exits to the T junction 
are long relative to this diameter; connecting lengths are approximately 100 mm. Flow rates for the incoming fluid are 
varied through the experiments, but typical rates are on the order of 1 ml/min which corresponds to fluid velocities 
in the T-junctions of 17 mm/s. The Reynolds number for water flow at these velocities is Re = 19. The Reynolds 
number of the sucrose solution of higher viscosity is reduced. The syringe pumps use glass syringes with a capacity of 
either 50 to 20 ml. All pumps were periodically tested and calibrated for their ability to deliver a steady and correct 
flow rate. 

In a typical experimental run, we wish to measure the volume fraction of fluid 2 in the two outlets as a function 
of the outlet flow fraction. The outlet flow fraction is defined as the flow in the branch divided by the total inlet 
flow; Qbr/Qin- We set the inlet pumps to a constant flow rate and fixed ratio to control the inlet volume fraction, 
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FIG. 1. Top view schematic of the experimental setup. Two syringe pumps push the two fluids at a controlled rate. Water is 
in inlet pump 1 and the viscous sucrose solution is in inlet pump 2. The volume fraction (of fluid 2) into the test T-junction is 
given as <3>i n = Q2HQ1 + Q2)- Fluid is collected from the branch and the volume fraction in the outlet is measured. The outlet 
pump can be switched to the branch and the experiment repeated. 



"Sin = Q^IQva-, and total flow rate Q ln . We then vary the outlet withdrawal flow rate of the pump attached to the 
run. After allowing sufficient time for the system to reach steady state, we collect the fluid from the branch, typically 
collecting at least 2 ml of fluid. The outlet solution is then well mixed in the collecting container and two successive 
samples are pipetted and measured with the Brixometer; all readings of the same collection sample fell within the 
stated error of the device (0.2 % error in mass fraction). We then change the flow rate of the outlet pump to collect 
data over a range of outlet flow fractions Qbr/Qin- After each change in the flow rate, we let the system come to 
equilibrium by waiting several flow through times for the network. 

With the outlet pump placed on the run, once Qbr/Qin becomes less than about 0.1 we cannot collect a sufficient 
amount of fluid from the branch before the input syringes empty, thus we end the experiment. We move the pump 
from the run to the branch and repeat the procedure, collecting and measuring fluid from the run while controlling 
the flow in the branch. In this configuration we can acquire data when the flow in the branch is small, but have 
difficulty obtaining data when Qbr/Qin > 0.9 since too little fluid accumulates from the run before the inlet syringes 
empty. 

While we could collect the contents of the branch (run) and determine the contents of the run (branch) from 
conservation, we avoid this approach. We measure the contents of the branch and the run in independent experiments 
and then confirm that the inferred values from conservation are in agreement with the measured values. It should be 
noted that the inferred values of the branch contents (when the run is measured) are unreliable when Qbr/Qin > 0.9. 
The sensitivity in the calculation is such that a small error in the measured value leads to a large error in the inferred 
value, thus the measured values arc to be trusted over the inferred. We repeat each measured data point three times. 
Each of these three measurements is a unique experiment; new inlet fluids are mixed, all syringes are cleaned, and a 
new identical network is constructed from new materials and components. 



B. Simulation 



We conducted simulations of the system using a commercial 3D finite element code, Comsol Multiphysics. The 
simulations solve the steady state, incompressible Navier-Stokes equations for conservation of momentum, 

pu- Vu= -VP + pg + V • (/i (Vu + Vu T )) , (1) 

and mass 

V • u = 0. (2) 

We solve the convection-diffusion of a dilute species which represents the relative concentration of sucrose, 

uVc = L>V 2 c. (3) 

Here, u is the velocity vector, p is the density, p is the viscosity, c is the concentration, and D is the diffusivity. The 
concentration equation couples back to the momentum equation through the dependence of the density and viscosity 
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FIG. 2. a) Simulation geometry. The color represents the concentration field. The two fluids are brought in vertically with 
respect to gravity so that the flow will be pre-stratified in the inlet, b) Sample velocity profiles for stratified flow as a function 
of radial distance from the center of the tube. Profiles are shown for viscosity ratios of 5, 15, and 100. For each case the inlet 
volume fraction is <I?i n = 0.5. Notice that the location of the interface between the two fluids moves toward the wall as the 
viscosity ratio increases. The radial coordinate in the plot is taken along a vertical line through the center of the tube in the 
direction of gravity. For a viscosity ratio of 1 we would have the classic parabolic velocity profile. 



M = M-) , (4) 



on concentration. We approximate the fluid viscosity as, 

i (c 2 +c)/2 

.Mi 

where p\ is the viscosity of water (fluid 1) and p2 is the viscosity of the inlet sucrose solution (fluid 2). The power 
of (c 2 + c)/2 was empirically fit to match the experimental data with reasonable accuracy & Our definition of 
concentration, c, is normalized by the concentration in fluid 2, thus the concentration in our model varies between 
and 1. The density is assumed to depend linearly on concentration, p = p\ + (p2 — Pi)c, where p\ and P2 are the 
densities of fluid 1 and 2 respectively. 

The simulation domain, shown in Figure UK, consists of two inlet tubes which are aligned with respect to gravity. 
The inlet concentration is in the upper inlet (pump 1) and 1 in the lower inlet (pump 2); both inlet flow rates are 
controlled. The two fluids meet and merge in a single tube, approach the test T junction and then proceed to the 
outlets. The stratified flow in the approach to the test T is enforced by aligning the inlet tubes with gravity. In the 
experiment, the long entrance length ensures fully stratified flow. The outlet flow rate is controlled on the branch 
through the boundary condition while the outlet at the run is allowed to freely flow out to a reservoir. The simulation 
was run on progressively finer grids to ensure adequate convergence in the solution. 

We make the equations dimensionless by using the pipe diameter, d, as the length scale, and the total input flow 
rate divided by the area of the pipe as the velocity scale, J7 . Under this scaling the equations become, 

(1 + /3c)u • Vu = -VP + Fr/3cz + ^-V • (/2 (c2+c)/2 (Vu + Vu T )j , (5) 

V • u = 0, (6) 

There are six dimensionless parameters: the Reynolds number Re = Pl j^° d ; the Froude number Fr = gd/U§; the 
viscosity ratio p, — p,\\ the inlet volume fraction $i„; the Schmidt number Sc = pj pD; and the density difference 
ratio (3 — {p2 — p\)j p\. Note that the Schmidt number is a material constant with a typical value of Sc = 1000 and 
is thus not a free parameter that we can easily control. 

When two fluids of different viscosity flow inside a tube, it is important to realize that the volume taken up inside 
the tube is not equal to the ratio of the flow rates. The low viscosity water is squeezed to the top of the tube where it 
has a higher velocity than the viscous sucrose solution. Sample velocity profiles in the tubes are shown in Figure [^p. 
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C. 



Results 



In Figure [3] we show a comparison of the experimental measurements and 3D simulations at four different total 
inlet flow rates. In each sub-figure we plot the volume fraction of the fluid in the branch and in the run, divided by 
the inlet volume fraction as we vary the flow ratio, Qbr/Qin- For all cases the inlet volume fraction is held fixed at 
<f>i n = 0.5 and ji = 15, but the overall inlet flow rate is varied. If there were no phase separation effect, then each 
plot would show unity for all values of the flow fraction. However, we see significant phase separation which depends 
strongly on the overall flow rate or inlet Reynolds number. Each point represents the average of three independent 
trials and the error bars represent the spread in these trails. The triangles represent the inferred values from the 
measurements of the other outlet. Across the range of experiments the measured and inferred values are in excellent 
agreement, giving us high confidence in the experimental procedure. 

In general, we find good agreement between the simulation and the experiment. The agreement is best at low flow 
rate while at higher flow the simulation under-predicts the experimental data. However, the predicted trend is in 
excellent agreement for a non-linear flow problem. Since the amount of separation depends sensitively on the total 
flow rate, it is clear that the non-linearity of the Navier-Stokes equations is the origin of the phase separation. The 
mechanism for the splitting is an inertia effect — as the fluid mixture approaches the T, the fast moving water tends 
to overshoot the 90 degree bend and continues onward along the run. The run always contains more water than the 
inlet fluid. 

The data near the end points (when Qbr/Qin is close to zero or one) are difficult to resolve both for the simulation 
and the experiment. For example, when the flow in the branch approaches zero, the volume fraction in the run must 
approach that of the inlet. The branch can have essentially any volume fraction and leave the run unchanged. In the 
experiment, we cannot collect enough fluid at these conditions before the syringe pump empties, thus the contents 
cannot be measured directly. We also cannot calculate the contents of the branch from the contents of the run, unless 
we know the run contents with extreme precision — the error gets amplified in this regime. Computationally, the flow 
and flux of concentration in the branch are both very close to zero near this point. Therefore, the ratio of the two is 
susceptible to a small amount of error in the overall numerical solution. With the current arrangement, we cannot 
confidently say what the concentration in the branch (run) is as the branch flow rate approaches zero (one). The 
simulation is well behaved in the range 0.05 < Qbr/Qin < 0.95. Beyond this range, while the numerical solution 
returns a well converged solution, the calculation of the volume fraction in the low flow outlet is not well behaved - 
exactly as in the experiments. 

In addition to varying the flow fraction, we also explored how the phase separation depends upon other parameters 
in the system. In Figure[l]we hold the outlet flow fraction fixed at 0.5 (i.e. Qbr = Qrun) and vary both the total flow 
and the viscosity contrast. The inlet volume fraction is $;„ = 0.5. In Figured we vary the overall flow rate while 
in Figure 3)d we vary the viscosity ratio. As before we measure the contents of the branch and the run separately, 
using three separate trials for each data point. The Comsol simulation tends to under- predict the phase separation, 
however the trend in both cases is well-captured. 

In Figure [5] we show the phase separation as we change the inlet volume fraction. In this case we hold the total 
flow at 4 ml/min and show results for $j n = 0.7 and <E>i n = 0.3. The viscosity ratio between the two fluids is /2 = 15. 
We find that the phase separation depends quite strongly on the incoming volume fraction, and that the separation 
is quite dramatic when the volume fraction of fluid 2 is low. 

In general we see find very good agreement between the experiment and simulations. Typically, the simulation 
under-predicts the amount of phase separation, however all trends seem to be well captured by the finite element 
simulation. The simulation therefore, is a useful tool for these flows. We can use the simulation to predict phase 
separation functions over a much wider range of parameter space, much more quickly than we can experimentally. 
Given the good agreement we find here, we have confidence that predicted separation functions will be accurate. To 
the best of our knowledge, these measurements and simulations represent the first in this geometry with laminar, 
stratified, miscible flow. 

As we will explain in the next section, these measured separation functions are a key ingredient in network models. 
The separation function must be applied at each node within the network. In order to quickly explore ramifications of 
phase separation for network flows, we empirically develop a simple one-parameter model function for <&b r and $ run . 
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FIG. 3. Comparison of experiment and simulation for the phase splitting at a T junction as a function of the overall flow 
rate. Each inlet pump has the same flow rate for the sucrose solution and the water which were a) 0.5 ml/min, b) 1 ml/min, 
c) 2 ml/min and d) 5 ml/min. Each data point represents the mean of six measurements (2 independent measurements on 3 
independent but identical networks). The error bars represent the maximum and minimum values measured. The composition 
of the branch and the run are measured separately, though this measurement is redundant. The triangles represent the mean 
value in the branch (run) as calculated from the measurements in the run (branch) — the data are consistent. The red data 
points are data measured in the run and the blue points are measured in the branch. The solid line is the finite element 
simulation from the branch and the dashed line is the simulation from the run. The sucrose solution coming from pump 2 is 
a 50% mass fraction solution which has a viscosity of ft — 15. The Reynolds number in the T-junction (based on the water 
viscosity) for each case is approximately a) 19, b) 38, c) 76, and d) 190. 



If the control parameter a is greater than 1 we use a piecewise approximation where $ run = in the regions where 
the volume fraction would be negative. If $ ran = is set to zero, then $br is computed from conservation. The 
construction of this simple function, while far from precise captures the basic trends of our experimental data. This 
empirical function allows us to scan parameter space easily in making network predictions. With this simple function, 
a is related to the strength of the phase separation. An example of the fit function compared to experimental data 
is shown in Figure [51 



III. NETWORK FLOWS 



Now that the phase distribution functions are known for a single junction, we can begin to use these functions to 
predict how the phases would distribute within a network. In this paper we use one of the simplest networks possible, 
Figure[?l which we demonstrated in previous work had bistability when the two inlets {Q ln ,i and Qi n ,i) were different 
fluids and the viscosity contrast was sufficiently strong (l6j . In this work we consider the same network, only the two 
inlets to the network are the stratified 2-component flow rather than single phase fluids. We find a different type of 
bistability which originates from the phase distribution function at the single T-junction, remarkably occurring when 
the two inlets to the network are identical fluids. In this setup, we use four inlet pumps to create a controlled stratified 
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FIG. 4. Comparison of experiment and simulation. Here the two outlet flows are fixed to always be equal; Qbr/Qin = 0.5 In 
a) we only vary the total flow rate holding the viscosity ratio at jl = 15. In b) we hold the total flow constant at 4 ml/min, 
Re = 76, and vary the viscosity ratio. The symbols are the same as in Figure [3] 




FIG. 5. Comparison of experiment and simulation for the phase splitting at a T junction. In both experiments the total inlet 
flow rate is 4 ml/min a) the sucrose inlet is 2.8 and the water is 1.2 ml/min, $i„ = 0.7 b) the sucrose is 1.2 and the water is 2.8 
ml/min, <E>i n = 0.3. The sucrose solution in the inlet pump has a viscosity 15 times that of water. The symbols are the same 
as in Figure |31 



flow as the two inlets to the network. The inlet volume fractions to the network are defined by $i n ,i = Q2HQ1 + Q2) 
and $i n! 2 = Qa/ (Q3 + Qi)- The total flow into the network is Qtot = Qin,i + Qin,2- 



A. Analysis 



The pressure drop across any length of tube can be described by 

AP = QR 

where R is the hydraulic resistance of the tube, AP is the pressure drop and Q is the volumetric flow rate. The 
resistance for a single phase laminar flow is simply Pouiselle's Law, R = 128/iL/7rei 4 where L and d are the length and 
diameter of the tube respectively. For our two-fluid system the viscosity in Pouiselle's law is replaced by an effective 
viscosity. The effective viscosity is calculated from the full velocity profile of the two fluids in contact with each other 
in the tube, as in Figure [2]}. In the flow regime of our network, the resistance behavior is similar to what one would 
calculate for immiscible fluids in fully developed flow [34|. The setup for the calculation of the effective viscosity of 
immiscible fluids is shown in Figure [5£i. Since the flow is fully developed along the axial length of the tube, we solve 
for the velocity across a cross section, of tube assuming a horizontal interface separating the two fluids. Since the 
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FIG. 6. Comparison of experimental data for the phase splitting at a T junction compared to the simple 1 parameter model 
over a range of conditions. The points are experimental data (error bars removed for clarity) and the lines is from the fit 
function - solid lines for the branch and dashed lines for the run. All experimental data is the same as shown in previous 
figures. The conditions shown are; i) black lines for a — 0.6 and black triangles for 2 ml/min total flow with $i n = 0.5; ii) red 
lines for a = 1.1 and red pluses for 4 ml/min total flow with <!>i n = 0.5; iii) blue lines for a — 1.5 and blue circles for 10 ml/min 
total flow with $i n = 0.5; iv) magenta lines for a — 2.2 and magenta diamonds for 4 ml/min total flow with <I>i n = 0.3; In all 
cases the sucrose solution in the inlet pump has a viscosity 15 times that of water. 
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FIG. 7. Top view schematic of the network as constructed in experiment. The results are sensitive to whether the connector 
between the two outlets is the run (as configured here) or the branch of the T-junction. 



velocity only varies across the cross section the Navier-Stokes equations simplify dramatically as illustrated in the 
schematic. Once the velocity profile is computed for a fixed interface location, we can easily calculate the overall 
hydraulic resistance and the volume fraction of the two flows. 

Since our system has miscible fluids, we can solve the same problem only for an interface which was been smeared 
by diffusion. When diffusion exists in the system, the effective resistance of the flow in the tube depends on how 
long the two fluids have been in contact with each other in the tube. For long and narrow tubes, there would be 
sufficient time for molecular diffusion and the behavior limits to the fully mixed case. For short tubes, we would limit 
to the immiscible result. The relevant dimensionless parameter would compare the time scale for diffusion across the 
cross section of tube relative to the transit time through the tube; i.e. ud 2 /DL — ReScd/L. Example curves for 
the effective viscosity (based on miscible fluids) are shown in Figure [5J these fall close to the immiscible result. The 
calculation of the effective viscosity and the effect of diffusion on the resistance in a single tube was discussed in more 
detail in our previous work[l6j. 

Once the effective viscosity as a function of the volume fraction is known, we can relate the pressure drop and flow 
in each of the tubes. For our network, the flow is defined as positive for flowing out at exits A and B and positive 
flow in C is defined from left to right. With this convention the pressure drops in the three tubes must satisfy, 

AP A = AP C + AP B 

Using flow conservation, we know that Qa = Qin.i — Qc and Qb = Qin.2 + Qc- We can describe the flow in tube C 
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FIG. 8. a) Schematic and formulation for calculating the effective viscosity of immiscible stratified flow in a tube. We show 
a cross section of the tube where u is the axial velocity coming out of the page. The boundary conditions are no-slip at the 
wall and constant stress at the interface between the two fluids, b) Effective viscosity for miscible stratified laminar flow of two 
fluids in a tube as a function of volume fraction of the more viscous fluid in the mixture. Sample curves are shown for viscosity 
contrast of 5, 10, and 15. 



as, 

Qin.lRA — Qin,2-RB 



Qc 



Ra + Rb + Rc 



It is important to remember that the resistances of each tube in the above expression depends upon the volume 
fraction in the tube. 

When the flow in C is positive then phase separation occurs at the T-junction on the left and thus this phase 
separation determines the contents of tubes A and C. The volume fraction in the three tubes follows, 

$A = ^brCTp^-jBfim,!, ( 8 ) 

*C = *ru„(7^-,Re in ,i,$ in ,i,/ii n ,i) (9) 
t<rin,l 

- Qi„,2$i„,2 + Qc$c 

®B = p: ■ (10) 

Qb 

When the flow in C is negative, 

$B = *br(^-,Rein.2,$in,2,Ain,2), (H) 
Win,2 

$C = $run(7^,Re ini2 ,$i„.2,Ain,2), ( 12 ) 
Win,2 

, Qin,l$in,l - QC$C ri „s 
®A = t: ■ (Id) 

Qa 

The functions $br and $ run are the phase separation functions which come from either experiment, 3D simulation, 
or our empirical one-parameter model. These functions, as described in Section II, depend upon the relative flow of 
the branch to the inlet, the Reynolds number, the inlet volume fraction, and the viscosity contrast between the two 
fluids. Once the volume fraction in the tubes is known, the effective viscosity of the network tubes are known from 
Figure [3] Given that the phase separation functions and resistance are not simple analytical results, a simple result 
cannot be derived. The problem is sufficiently defined to numerically find the volume fraction and flows through the 
entire network for fixed inlet conditions. The above equations provide the non-linear relationship for Qc/Qtat as a 
function of Q in ,i/Qtot- 

Some sample predictions using the one-parameter model are shown in Figure [S] for an inlet volume fraction of 
$in.i = ^in.2 = $in = 0.5 and a network where all tubes are equal lengths. In Figure we show the effect of 
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FIG. 9. a) Predictions of the flow in tube C as a function of the ratio of flow in inlet 1 to the total flow rate Qm,i/ (Qm,i + Qm,2) 
for values of a = 0, 0.5, 1, and 2; (blue, red, black, and magenta respectively). The viscosity contrast is fixed at fi = 150 
b) Predictions of the flow in tube C for different viscosity contrast of jl = 1, 4, 20, and 150; (blue, red, black, and magenta 
respectively). The parameter for the phase separation function is a = 2. In both figures a) and b), the inlet volume fraction is 
"I'm,! = ^in^ = 0.5 and the lengths of all tubes in the network are the same. 

variation in the one parameter in the phase separation function, a, at a high viscosity contrast between the two fluids, 
jl = 150. It is clear that if sufficiently strong phase separation occurs then bistability in this simple network can 
emerge. For the curves a — 1 and a = 2 there are three possible values of Qc when Qi nj i/Qtot = 0.5. The value at 
Qc — is not stable and thus would not be observed in experiments. While bistability exists for a = 1, it is much 
more extreme for larger values of a. In Figure [Hb, we show the effect of the viscosity contrast between the two fluids 
in the system for a fixed value of a = 2. We see that the window of bistability grows with the viscosity contrast. 
While these figures do not span all possible parameters, it is clear that for bistability to emerge in these networks the 
phase separation and the viscosity contrast must both be sufficiently strong. 

We can easily derive a criterion for the emergence of bistability in this network flow. We can take the formulation 
for the network and analyze the derivative of the curve dQc I dQin,i around the equilibrium point when Qc = 0. At 
this trivial equilibrium point Qi n! i = Q; n ,2 and Qa = Qb- Looking at Figure [TUJ we see bistability emerge when the 
slope at this point becomes vertical. Considering the special case where the two inlet pumps to the network have 
identical fluids and inlet volume fractions (i.e. $j n) i = $i n .2 = Sin and Mi n ,i = Min,2 = Min)> we obtain the slope of the 
flow curve to be, 



8Qc 

OQin 



1 



1 



vc(Qc=Q) 



Lb+La 



,)■{' 



(*in - ®C(QC = 0)) 



(14) 



Here, nc and <&c are the viscosity and volume fraction in C when Qc — 0. The value of $c(Qc — 0) comes from 
the phase separation function. For the network configuration shown in Figure O tube C is the run of the T-junction. 
For simplicity of the expression, we have assumed that all tube diameters are equal in Equation 1141 

We can easily see that bistability would emerge when the denominator of Equation [14] equals zero. Without loss of 
generality, we assumed through our definitions that viscosity must increase with volume fraction and thus, ^ > 0. 
Since the tube lengths and viscosity must all be positive numbers, Equation Q3] tells us some general things about 
when bistability cannot exist. In the limit of no phase separation, $c — ^im and the denominator would always 
be positive and bistability could not exist. Likewise, if the phase separation is such that $c(Qc = 0) > <E>i n , then 
bistability is also not possible. If $c < Sj n then bistability is possible depending on the network geometry and 
viscosity of the fluids. The criterion for the onset of bistability can be stated as, 



Min 9$ 



yc(Qc=0) L c 



$i„-$c(Qc = 0) 



(15) 



To obtain better intuition about when this criteria is met, let us suppose that the effective viscosity follows a simple 
law, fi = where jl is the viscosity contrast between the two fluids in the system. With this viscosity law, 
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■jjj- = log(/i). Therefore the required viscosity contrast between the two fluids, /2, would be exponential in 

the right hand side of Equation 1151 While our viscosity law is different it has a similar functional form, thus a small 
change in the right hand side will get magnified by exponential behavior. 

In the limit of Lc > > La + Lb the right hand side of Equation [15] would grow and thus bistability would become 
very unlikely. In the limit when Lc << La + Lb the criteria would simplify to, 



_L 2t 
Min d<Z> 



1 



$in - MQc = o) 



This value would represent the minimum viscosity contrast needed to observe bistability in any arbitrary network. 

We can test Equation [15] against Figured where $i n = 0.5. In this figure La = Lb = Lc and the phase separation 
is strong such that <&c(Qc = 0) = 0. For these conditions and our stratified viscosity law, solving our criteria 
numerically (there is no analytical expression) would provide a critical viscosity contrast of fl = 26. In the Figure [Hb 
we show a curve for fl = 20, and the slope is nearly vertical, though if we zoom in we find that it is still in the positive 
direction. 



B. Experiment 

A schematic of the experimental network setup is shown in Figure [7] We use four inlet pumps to create a controlled 
stratified flow as the two inlets to the network. In order to easily monitor the flow rate in the exit branches of the 
network, we simply leave exits A and B open to atmospheric pressure. With this simple setup we collect the fluid 
from both outlets over a fixed time period in order to measure the exit flow rates, Qa and Qb- The exit flows are 
collected in a disposable test tube and the fluid mass is weighed on an analytical balance. 

In the network experiments, we used the same T-junction geometry described in Section II. However, in these 
experiments after the T-junction the tubes A, B, and C were narrowed down to tubing of 0.5 mm inner diameter. 
This decrease in diameter was done to increase the hydraulic resistance between the network connections. Since the 
outlets A and B were held open at atmospheric pressure, we want the hydraulic resistance of the tubes to significantly 
dominate over surface tension effects due to dripping at the exits or hydrostatic pressure differences due to A and B 
not being at precisely the same height. In all experiments reported here the hydraulic resistance of A, B, and C were 
the same. 

In these experiments the pumps on the left and right side are always kept in the same proportion such that the 
volume fraction in the two inlets are constant at $i n ,i = 3?i n ,2- In the experiment the pumps on inlet 1 in Figure [7] arc 
held at a constant flow rate of 2 ml/min. The total flow rate in inlet 2 is varied between 8 and 0.5 ml/min in order to 
change the inlet network flow ratio, Qi n ,i/Qtot- Each experimental data point is repeated in three independent trials 
where new networks are constructed and new solutions mixed. 

We used our model phase separation function to help guide us to regimes where strong bistability would be expected. 
We determined from our experiments that strong phase separation occurs when the inlet volume fraction is 30%, see 
Figure [B^,. We also determined from our network analysis that a viscosity contrast of 150 would give us a strong 
window of bistability in the network; Figure[9}D. Since we found in earlier sections that the amount of phase separation 
is not very sensitive to the viscosity contrast, we assume that the fit value of a — 2.2 that was found in Figure^ 
for a viscosity contrast of 15 is approximately equivalent to that at a contrast of 150. We use the fit function with 
a = 2.2 in our network prediction. 



C. Results 



The predictions of the network model compared to our experiments are shown in Figure [TTJk . We show the prediction 
using the simple phase separation fit function; shown as the solid line. To further validate our predictions, we also used 
our finite element simulation at the experimental parameter values and get the predicted phase separation function 
as well (the result is very similar to Figure [SJa). The network prediction based on the phase separation function 
calculated using the 3D finite element code is shown as the dashed line. A clear window of bistability exists and 
the behavior is quite remarkable. If we slowly adjust the incoming flow rates around the point where Qi n ,i — Qin,2 
then the flow in the outlet undergoes a dramatic jump when we pass the transition point. The agreement between 
experiment and theory are excellent. Note that the window of bistability predicted by the model using the phase 
separation which came from the finite element simulation is not as strong as we find experimentally. This result is 
consistent with what we found in an earlier section where the finite element simulation typically under-predicts the 
amount of phase separation seen in experiment. 
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FIG. 10. a) Flow in exit branch A of the network in Figure [TJd determined from experiment and theory. Notice the region of 
bistability around the state where the flow rates on the inlets are equal. Both inlets are a 30% sucrose and 70% water mixture 
in stratified flow. The viscosity contrast between the sucrose solution in inlet pumps and water is 150. The inset shows the 
network topology. The error bars on the red experimental data points show the maximum and minimum value recorded in 
three independent trials and the data points shows the mean. The blue data points (and error bars) were taken on two of the 
trials where we proceeded from left to right along the flow curve. The black data points were from one trial where we moved 
from right to left. In b) we have the exact same setup, only the connection for the network is made by the branches of the 
two inlet T-junctions - see the inset. In both figures the points are experimental data, the solid line is from the theory, using 
the fit function with a = 2.2 and the dashed line is from the theory using the phase separation as predicted from the Comsol 
simulation. 



In the experiment the pumps on inlet 1 in Figure [7] are held at a constant flow rate of 2 ml/min. When the system 
is in a state such that the flow in C is positive (generally Qi n ,i/Qtot > 0.5) then phase separation only occurs at 
the T-junction on inlet 1. Since this flow rate is constant, we expect a single value of a to suffice for describing 
the behavior. When the flow is in a state where the flow in branch C is negative (generally Qi n ,i/Qtot < 0.5) then 
phase separation only occurs at the T-junction between branches B and C. In order to span values of Qi n ,i/<9tot, the 
pumps on Inlet 2 are varied from an inlet flow rate of 8 to 0.5 ml/min in the experiment. Since the phase separation 
function is sensitive to the total flow rate, then the value of a would be expected to be variable. The experimental 
data show an asymmetry about <5i n ,i/Qtot = 0.5 which is not present in the simple model which assumes constant 
phase separation behavior. Despite these differences, the agreement between the data and the prediction is excellent. 
Note the phase separation function calculated with the 3D finite element simulation was run at a constant total flow 
rate of 4 ml/min to approximately correspond to the average total flow in the network across the range of parameters. 

Another interesting feature of this system is that if we keep the same basic network but reconfigure the system 
such that tube C in the network is the branch of both inlet T-junctions rather than the run, the window of bistability 
disappears completely. These data and the predictions are shown in Figure [TOb . The experimental data are consistent 
with the model. Again the data show an asymmetry about Qi n ,i/Qtot which not present in the model due to our 
assumption of constant phase separation for all inlet flow rates. 

An interesting feature of Figure [TUb. is that both the data and the prediction using the phase separation from the 
finite element simulation show a flattening of the slope around the equilibrium point, Qc = 0, that is not seen with 
the simple fit function. The reason is that the fit function assumes that the volume fraction in the branch goes to 
$i n when the flow in the branch is zero; ie. 3>br(Qbr = 0) = $i n = 0.3. This assumption is made to keep the fit 
function to one adjustable parameter, the assumption is not based on data or some physical mechanism. As explained 
in Section II, the experimental setup does not allow us to accurately determine the branch volume fraction at this 
point. However, as can be seen in the phase separation functions computed by the finite element simulations show 
that the volume fraction in the branch as the flow in the branch goes to zero is likely somewhat less that the inlet 
volume fraction. For the case here, the finite element simulation predicts 3>br(Qbr = 0) ~ 0.18. Since the flow in tube 
C is close to zero when Q ln ,i/Qtot ~ 0.5, it is the phase separation behavior when the branch flow is close to zero 
which is important. 

What is interesting about the dashed curves in Figure [TU] is that the prediction is based entirely on first principles 
and simulation. The complex phase separation functions and effective viscosity laws are calculated with a Navier 
Stokes simulation. The dashed curves in these figures have no fit parameters or assumptions inherent in them. The 
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solid curve is based on a simple one parameter phase separation model that is convenient for exploring behavior and 
parameter space, but it should not be taken as more than a fit function. 

Finally, the criterion in Equation [15] provides a simple explanation as to why bistability is seen in one network 
configuration but not the other. In the configuration of the network shown in Figure 110a . tube C is the run of inlet 
T-junctions and thus $c(Qc = 0) « 0. Under these conditions bistability is possible when ft > 50; our viscosity 
contrast of 150 exceeds the minimum value. For the configuration where the network is connected such that tube C 
is the branch of the T-junction, $c(Qc = 0) = if we use the fit function; bistability is not possible under this 
condition. If we use the value of $c{Qc = 0) = 0.18 which is computed in the finite element simulation, then the 
required viscosity contrast would be jl > 5900. 



IV. CONCLUSIONS 

In this work, we have demonstrated that for laminar flow of two miscible fluids of different viscosity approaching 
a T-junction has significant phase separation. We measured phase distribution functions for this system which, to 
the best of our knowledge, has not been previously done. We compared our experiments on these phase distribution 
functions to 3D Navier Stokes simulations and find excellent agreement. We find that in the geometry used in this 
work, that inertia was the mechanism behind the phase separation. We described the phase separation as a function 
of the key parameters such as overall flow rate and viscosity contrast between the two fluids. 

Once such phase separation functions are known, we began to explore what the consequences would be for phase 
distribution within a simple fluid network. We find theoretically and experimentally that the unequal separation of 
two phases at a single network bifurcation can lead to multiple stable equilibria in a simple network. Even in the 
simple network studied in this work, we find that multiple equilibrium states can be observed. We derive a simple 
criterion for the existence of multiple equilibrium states and found to this criterion to explain the experimentally 
observed behavior. The phase separation functions depend on the exact details of the system; the fluids used and the 
geometry of the network bifurcation. However, the network results are generic and could be applied to or found in 
different systems. 
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